#include "KnotVectors.h"

namespace cagd
{
    // special constructor
    KnotVector::KnotVector(Type type, GLuint p, GLuint n, GLdouble u_min, GLdouble u_max):
            _type(type),
            _order(p),
            _control_point_count((type == PERIODIC) ? n + p : n + 1),
            _u_min(u_min), _u_max(u_max)
    {
        switch (_type)
        {
        case UNCLAMPED:
        case PERIODIC:
            {
                _knots.ResizeColumns(_control_point_count + _order);

                GLuint   m    = _control_point_count - 1;
                GLdouble step = (_u_max - _u_min) / (m + 2 - _order);

                for (GLuint j = 0; j <= m + _order; ++j)
                {
                    _knots[j] = _u_min + (j - (GLdouble)_order + 1) * step;
                }

                _knots[_order - 1] = _u_min;
                _knots[m + 1] = _u_max;
            }
            break;

        case CLAMPED:
            {
                _knots.ResizeColumns(_control_point_count + _order);

                for (GLuint j = 0; j <= _order - 1; ++j)
                {
                    _knots[j] = _u_min;
                }

                for (GLuint j = n + 1; j <= n + _order; ++j)
                {
                    _knots[j] = _u_max;
                }

                GLdouble step = (_u_max - _u_min) / (n + 2 - _order);

                for (GLuint j = _order; j <= n; ++j)
                {
                    _knots[j] = _u_min + (j - (GLdouble)_order + 1) * step;
                }
            }
            break;
        }
    }

    // get knot value by value
    GLdouble KnotVector::operator [](GLuint index) const
    {
        return _knots[index];
    }

    // get knot value by reference
    GLdouble& KnotVector::operator [](GLuint index)
    {
        return _knots[index];
    }

    // determines an index stored in variable i for which u is in [u_{i}, u_{i + 1})
    GLboolean KnotVector::FindSpan(GLdouble u, GLuint& i) const
    {
        if (u < _u_min || u > _u_max || _control_point_count <= 1 || _knots.GetColumnCount() < _control_point_count + _order)
        {
            return GL_FALSE;
        }

        if (u == _knots[_control_point_count])
        {
            i = _control_point_count - 1;
            return GL_TRUE;
        }

        GLuint left  = _order - 1;
        GLuint right = _control_point_count;

        i = (left + right) / 2;

        while (!(u >= _knots[i] && u < _knots[i + 1]))
        {
            if (u < _knots[i])
                right = i;

            if (u >= _knots[i + 1])
                left = i;

            i = (left + right) / 2;
        }

        return GL_TRUE;
    }

    // evaluates non-vanishing normalized B-spline function values in a lower triangular matrix
    GLboolean KnotVector::EvaluateNonZeroBSplineFunctions(GLdouble u, GLuint &i, TriangularMatrix<GLdouble>& N) const
    {
        if (!FindSpan(u, i))
        {
            return GL_FALSE;
        }

        N.ResizeRows(_order);

        N(0, 0) = 1.0; // N_{i}^{1}(u)

        GLdouble left, right;

        // N_{j}^{k}(u), j = i-k+1, i-k+2, ..., i
        for (GLuint k = 2; k <= _order; ++k)
        {
            // j = i - k + 1
            right = _knots[i + 1] - _knots[i - k + 2];
            N(k - 1, 0) = (right > 0.0) ? N(k - 2, 0) * (_knots[i + 1] - u) / right : 0.0;

            // j = i
            left = _knots[i + k - 1] - _knots[i];
            N(k - 1, k - 1) = (left > 0.0) ? N(k - 2, k - 2) * (u - _knots[i]) / left : 0.0;

            // j = i - k + 2, ..., i - 1
            GLuint offset = i - k + 1;

            for (GLuint j = i - k + 2; j <= i - 1; j++)
            {
                left  = _knots[j + k - 1] - _knots[j];
                right = _knots[j + k] - _knots[j + 1];

                N(k - 1, j - offset) =
                        ((left  > 0.0) ? N(k - 2, j - offset - 1) * (u - _knots[j]) / left : 0.0) +
                        ((right > 0.0) ? N(k - 2, j - offset) * (_knots[j + k] - u) / right : 0.0);
            }
        }

        return GL_TRUE;
    }

    // get order of the knot vector
    GLuint KnotVector::GetOrder() const
    {
        return _order;
    }

    // get type of the knot vector
    KnotVector::Type KnotVector::GetType() const
    {
        return _type;
    }

    // get control point count
    GLuint KnotVector::GetControlPointCount() const
    {
        return _control_point_count;
    }
}
